(N 
> 



Asynchronism Induces Second Order Phase 
Transitions in Elementary Cellular Automata 



Nazim Fates 



O 
O 

(N 

■ "Nazim . FatesQloria . f r 

[Jh : LORIA - INRIA Nancy Grand-Est 

^ \ Campus Scientifique B.P. 239 

54 506 Vandoeuvre-les-Nancy, France. 

February 14, 2008 



Abstract 



Cellular automata are widely used to model natural or artificial sys- 
tems. Classically they are run with perfect synchrony, i.e., the local rule 
is applied to each cell at each time step. A possible modification of the 
updating scheme consists in applying the rule with a fixed probability, 
called the synchrony rate. For some particular rules, varying the syn- 
chrony rate continuously produces a qualitative change in the behaviour 
, of the cellular automaton. We investigate the nature of this change of be- 

' haviour using Monte-Carlo simulations. We show that this phenomenon 

^ is a second-order phase transition, which we characterise more specifically 

. as belonging to the directed percolation or to the parity conservation uni- 

' versality classes studied in statistical physics. 



keywords: asynchronous cellular automata, stochastic process, discrete dy- 
namical systems, directed percolation, parity conservation, phase transitions, 
■ universality class, power laws 

X 

\—i ' Foreword: In this article, we only present a limited set of plots showing our 

numerical experimentation. The complete set of graphs can be accessed at: 
|http: //www. loria. f r/"^f ates/Percolation/results .html 

The simulations were obtained with the FiatLux CA simulator |9|. 
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1 Introduction 



With the increase of computing power, ceUular automata (CA) are becoming a 
popular tool used to simulate various real-world systems. While early research 
was mainly concerned with the study of logical properties of abstract models 
|24| , many efforts now focus on using CA which closely mimic natural or artificial 
phenomena. 

Our research aims at studying the robustness of cellular automata to asyn- 
chronous updating, i.e., at evaluating to which extent a small modification of 
their updating scheme may perturb their behaviour. To tackle this problem, 
we propose to study not only a single model but a family of models, where the 
members are obtained by varying the updating scheme PU] and keeping the 
local rule constant. One simple way of producing such variations is to consider 
the so called a- asynchronous dynamics [14j . in which each cell updates its state 
with probability a, the synchrony rate, independently at each time step. 

To our knowledge, the problem of comparing the behaviour of synchronous 
versus asynchronous CA was first addressed by means of simulation in [6^ , with 
a qualitative evaluation of the changes. Other experimental studies followed 
showing that the updating scheme was indeed a key point to study [3l |3T1 [29] . 
On the theoretical side, few results have been obtained so far: the independence 
on the order of update history was studied in [151 ES] , existence of stationary 
distributions for infinite systems was studied in |21j and a first classification 
based on the convergence time was proposed in [131 114j . 

The Qf-asynchronous updating was studied experimentally [12] and it was 
shown that the 256 Elementary Cellular Automata (ECA) produced various 
qualitative responses to asynchronism. We were surprised to observe that, in 
some cases, a small variation of the synchrony rate a could trigger a qualita- 
tive change of behaviour. The present work is devoted to understanding these 
transitions with an experimental approach. We extend the first investigations 
presented in [llj . examining a larger set of rules with an improved precision. 

In |12j . we conjectured that the origin of the qualitative change of behaviour 
was a phase transition that belonged to the universality class of directed per- 
colation (see below). We emphasise that this hypothesis was mainly supported 
by the visual observation of the space-time diagrams produced near the critical 
point. The purpose of this article is to investigate this hypothesis numerically 
with large Monte-Carlo simulations. 

2 First observations 

This section introduces formal notations and presents a first set of simple ex- 
periments. 
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2.1 Formal definitions of the model 



Let a ring of n cells be indexed by £ = TLjnL ; a configuration is an assignment 
of state, or 1, to each element of £ ; the space of configm-ations is {0, 1}'''. The 
density of a configuration is the ratio of cells in state 1 over the configuration 
size n. The kinks density is the ratio of the number of 01 or 10 patterns over 
the size of the configuration n. An Elementary Cellular Automaton (EGA) is 
described by a function / : {0, 1} {0, 1} called the local rule. EGA are 
indexed according to Wolfram's usual notation [32 • 

We define the a- asynchronous updating scheme as the operation that consists 
in applying the local rule / with a probability a or keeping the same state with 
a probability 1 — a, independently for each cell. This updating scheme defines 
a probabilistic global rule which operates on the random configurations {x*)ten 
according to: 

1 ^li ^l+i) with probability a 
' * 1 x* with probability 1 — a 

By taking a = 1, we fall back on the classical synchronous case and as a is 
decreased, the update rule becomes more asynchronous while the effect of an 
update remains unchanged. For < a < 1, x* is a random configuration that 
depends on the sequence of cells that are updated at each time step. 



2.2 Selecting the Rules to Study 

In our first work [T^] , we experimentally detected a first set of rules that showed 
a qualitative change of behaviour when a was varied continuously. We re- 
examined this change of behaviour for the 88 Minimal Representative EGA in 
more details. For each rule, we arbitrarily fixed the ring size to n = 10 000 and 
we varied a with an increment of 1% from 0.02 to 1. For each a, we started 
from a uniform random initial configuration and we measured the evolution of 
the density during IGOGG/q; steps, putting a limit of 100 000 steps for a < 0.1. 
We extracted the average over the second half of the sample, the first half being 
used as a transient period to ensure that the system was stabilised. This average 
density can be considered as an estimate of the stationary density, i.e., the limit 
density that would be reached if the system size and the transient time grew to 
infinity. 

Figure [1] shows the estimated stationary density (or kinks density) versus a 
for EGA 6, 50, and 178 and Figure [2] shows how the variation of a affects the 
space-time diagrams these rules. We observe that for EGA 6 and 5G (respec- 
tively EGA 178), a non-zero density (resp. kinks density) corresponds to the 
existence of stable branching- annihilating patterns and a zero density (resp. 
kinks density) corresponds to a rapid extinction of these patterns. 

Exploring systematically the 88 Minimal Representative EGA, we observed 
that EGA 18, 26, 58, 1G6, and 146 all have plots which are similar to EGA 5G. 
Their behaviour is compatible with a second order phase transition: there exists 
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EGA 50 
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Figure 1: (top) Stationary density versus synchrony rate for EGA 6 (left) and 
EGA 50 (right), (bottom) EGA 178 : stationary density (left) and kinks density 
(right) versus synchrony rate. The discontinuity for EGA 6 and 178 at a = 1 is 
not an artifact. 

a critical value of a such that for a > ac, the stationary density is non-zero (the 
active phase) and for a < the stationary density is zero (the inactive phase) . 
The curve is continuous and reaches zero with an infinite slope for the critical 
value OLc- As our hypothesis is that the nature of the phase transition is directed 
'percolation (see Section [3.2p with an active phase obtained for high values of a, 
we call these rules the DPhi EGA. 

EGA 6, 38 and 134 displayed a similar behaviour but what is more surprising 
is that their phase transition is in an "inversed" pattern: the active phase is 
obtained for small values of a, the inactive phase for large values of a. We call 
these rules the DPiow EGA. 

The case of EGA 178 is also peculiar since the density curve is stable for large 
values of a but does not stabilise value for small values of a. This indicates that, 
for this particular rule, the choice of measuring density is not suitable. Instead, 
if we examine the variation of the kinks density, we obtain a smooth curve with 
a discontinuity compatible with a second order phase transition. We call this 
rule the DP2 EGA. The choice of this name is justified in the next section. 



4 



a = 0.25 



a = 0.75 



a = 0.75 



< 
O 
H 



a 



o 

o 





a 



00 

1> 



o 





a 




Figure 2: Space time diagrams for EGA 6 (top), EGA 50 (middle) and EGA 178 

(bottom). Synchrony rate is varied : a. = 0.25 (left), a = 0.50 (middle), a = 0.75 
(right). Time goes from bottom to top; the time factor is rescaled by a factor 
l/a(i.e., for a = 0.25 only time steps that are multiples of 4 are displayed). 
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3 Phase Transitions 



In this section, we present a short review of works related to asynchronous or 
probabiHstic ceUular automata and second order phase transitions. 

3.1 Universality classes 

Many physical or numerical systems exhibit critical phenomena: a continuous 
change in the value of a control parameter may produce a discontinuous response 
at the macroscopic level. It is remarkable that near the critical transition, 
the laws governing the evolution of such systems are generally power laws (see 
below). In short, one may understand the origin of these power laws from the 
fact that near the transition point, the system has a self-similar fractal structure 
and no typical spatial wavelength. An important property is that the same 
exponents of the power laws, the critical exponents, may be found for different 
models that do not necessarily share the same definitions at the microscopic 
level. 

The collection of all models that are described by the same set of critical 
exponents is called a universality class. Looking at the space-time diagrams 
produced by the asynchronous EGA, we originally found similarities with the 
patterns produced by couple map lattices used in hydrodynamics [T]. These 
similarities lead us to identify our phase transition as possibly belonging to the 
directed percolation (DP) universality class [12j . 

It is beyond the scope of this paper to list all the models that belong to the 
directed percolation universality class and we refer to [T7l[26] for a review. As far 
as cellular automata are concerned, directed percolation was observed in various 
problems. To our knowledge, the first CA that was shown to exhibit DP is the 
Domany-Kinzcl cellular automaton f20', '8] . This model is a tunable probabilistic 
CA that has the ability to display two different transitions depending on the 
tuning of its parameters. Various other models involving probabilistic CA were 
also shown to exhibit DP phenomena {e.g., p71 1^ ). 

Directed percolation was also identified in problems involving synchronisa- 
tion of two instances of cellular automata {e.g., [161 [30]). To our knowledge, 
the only example of directed percolation induced by asynchronism was given by 
Blok and Bergersen for the famous Game of Life 4 . The protocol they used to 
identify the universality class of the phase transition relied on the measure of a 
single critical exponent, the (3 exponent (see below). 

3.2 Directed Percolation 

Percolation problems are studied in the fields of discrete mathematics and sta- 
tistical physics. They were initially motivated by the need to model situations 
in which a fluid evolves into a porous random medium [5] . In the classical prob- 
lem of isotropic percolation, the porous medium is modelled by a regular infinite 
two-dimensional square lattice. The nodes of the lattice, or .sites, can be either 
open (with probability q) or closed (with probability I — q) ; the links between 
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Figure 3: Example of directed bond percolation: filled circles are wet sites, 
empty circles represent dry sites ; solid (resp. dashed) represent open (resp. 
closed) bounds. The links can be open with a probablity p. The problem is 
to determine the critical value Pc for which the size of the clusters of wet sites 
diverges. 

the sites, or bonds can also be either open (with probability p) or closed (with 
probability 1 — p). Starting from an initial set of wet sites, the question is to 
determine the set of sites that will also be wet if the liquid flows into open sites 
and open bonds. The sets of all wet sites connected through closed bonds and 
sites is called a cluster. 

If we set q = 1 (respectively p — 1), we have the isotropic bond (resp. site) 
percolation. For the isotropic bond percolation problem, the average cluster 
size diverges for the critical value Pc = 1/2 : we say that the fluids percolates 
through the medium. Directed bond percolation is a non-isotropic variant of the 
previous model in which links are oriented according to a particular direction 
(see Figure |3]). This models situations in which the fluid can go only in one 
direction, for example when submitted to gravity. 

We can also formulate directed percolation as a probabilistic dynamical sys- 
tem: in this case, time plays the role of the non-isotropic dimension. More 
formally, if we represent the state of a site i e N at time t by s* € {0, 1} (dry or 
wet site), starting from an initial condition , the states of sites are updated 
according to the simple rule [T7j : 




1 if [sti = 1 and Clip) = 1] or [sj^^ = 1 and 7^*(p) = 1] 
otherwise 



where (£*) and (JZl) are i.i.d. Bernoulli random variables ; they model the 
probability for the left or right bond to cell {i,t) to be open or closed. 

Theory and observations predict that for an infinite lattice and for a fixed 
value of p, if we start from an initial configuration with all sites in wet state, 
the density of wet sites d(p, t) evolves to a positive limit for p > pc and to a 
zero limit for p < Pc T7]. More precisely, if we denote by doo(p) the infinite 
time limit of d{p, t), for p > pc, near the critical point {p pc), the asymptotic 
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density do^ diverges from zero by following a power law: 

c^oo(p) (p-Pcf 

Note that as we have /? < 1, the right derivative of doo{p) has an infinite value 
at the critical point. This explains why the transition is experimentally seen as 
"abrupt" when p is varied smoothly. At the critical point p = Pc, the density 
vanishes to zero doo(Pc) = and the decrease follows a power law: 

d{p,,t)^r^ 

The values of the two critical exponents 5 — 0.1595 and /3 — 0.2765 are 
known by numerical simulations (the values are given here with four digits, see 
|17| for a better precision). Determining their values analytically is difficult 
and it is so far an open problem to know whether these exponents are rational 
numbers. There are also other critical exponents that can be used to identify a 
universality class but we choose here to focus only on the measure of /3 and 5 
exponents (following [5Sl[Tnj for instance). 

3.3 The PC-DPa class 

In the case where the local rule has two symmetric states, it is generally ob- 
served that the phase transition is either in the Z2 -symmetric directed percola- 
tion (DP2) class or in the parity- conserving class (PC). Again, we refer to pTll^S] 
for a detailed description of these two universality classes, their similarities and 
differences. In short, the PC class appears when one considers models with 
branching-annihilating random walks with an even number of offspring {e.g., 
|19|1. The DP2 class is observed with models that introduce two symmetric 
states, as it is the case for EGA 178. 

It is well known that the two classes coincide for dimension 1 and differ for 
higher dimensions. An intuitive reason for this surprising property is indicated 
by Hinrichsen [T7] : "active sites of PC models in c? > 2 dimensions can be 
considered as branching-annihilating walkers, whereas DP2 models describe the 
dynamics of branching annihilating interfaces between oppositely oriented in- 
active domains". From the observation of space-time diagrams of EGA 178 and 
by analogy with other models, we conjectured that the phase transition of this 
rule belonged to the DP2 universality class [TU] . 

3.4 Hypotheses to test 

The phase transition theory stipulates that there should be two macroscopic 
parameters, respectively called the control parameter and the order parameter, 
that satisfy the power laws near the critical point. For the sake of simplicity, 
we use the synchrony rate a as a control parameter and the density d as an 
order parameter. However, note that other possibilities may also be examined, 
for example in ^ the authors also use the activity [i.e., the ratio of cells in 
an unstable state) as an order parameter. For the nine DP EGA identified in 
Section [21 we thus expect to measure: 
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• d{ac,t) ^ t ^ at the critical point, 

• doo{a) ^ near the critical point, for the active phase, 

with Aq — a — Uc for the DPhi EGA and ~ ac — a for the DPiow EGA. For 
the DP2 rule EGA 178, we take the kinks density as an order parameter. By 
contrast with the DP universality class the exponents of the PC and DP2 classes 
are known analytically. In particular we have (5dp2 = 2/7 for two-dimensional 
lattices. (Recall that space and time play a role analogous to the two dimensions 
of the grid in the original directed percolation problem.) 



3.5 Protocol and Measures 

The measure of the DP critical exponents is a delicate operation that generally 
requires a large amount of computation time. The main difficulty resides in 
avoiding systematic errors when obtaining statistical data near the transition 
point. Authors have been mislead by their measures and concluded that a 
phase transition phenomenon was not in the DP universality class, which has 
later been proved wrong by using a different protocol and more precise measures 
[TSKinj. In order to limit the influence of systematic errors, we take the two-step 
protocol used in [16]: 

• We measure the critical synchrony rate ac by varying a until we reach 
the best approximation of a power law decay for the density. This first 
experiment also allows us to measure the critical exponent d. 



We measure the stationary density dao as a function of |a — ad and then 
fit a power law in order to calculate p. 



Note that these two steps are not independent since the second operation uses 
the previously computed value of ac- 

The other main concern is to evaluate the influence of finite size effects 
and metastability. In the active phase, finite DP systems are in an out-of- 
equilibrium state. This means that although infinite-size systems have a stable 
probability measure for the distribution of states, the finite-size systems used 
in the simulations may attain the absorbing state even if they are in the active 
phase. There exist many techniques to handle this problem, for example adding 
a small amount of noise to prevent the system from staying indefinitely in the 
absorbing state. In this work we choose to use large size lattices and verify 
experimentally that the results are not influenced by the trajectories that touch 
the absorbing state (here 0'''). 

We will present the curves for EGA 50 while only numerical data will be 
given for the other DP EGA. Indeed, rule 50 can be written in the rather simple 
form: 



V(a,6,c) G {0,l}^ f{a,b,c) 



if (a,6,c) ^ (0,0,0) 
if (a,6, c) = (0,0,0) 
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It is possible to see this rule as a one dimensional version of an "epidemic" rule: 
a healthy cell (state 0) gets infected (state 1) if at least one of its neighbour is 
infected ; once it is infected, it becomes healthy at the next update of the cell. 

4 Finding the Critical Synchrony Rate ac 

Following our protocol, we estimate the critical point by determining the change 
of convexity of the density curves in a log-log plot. The concave function char- 
acterises the active phase as the asymptotic density evolves towards a non-zero 
value ; the convex function corresponds to the inactive phase as the density 
evolves to zero with an exponential decay. 

For EGA 50, the previous experiment (Figure [1]) allowed us to locate the 
critical synchrony rate around ac ~ 0.62. Figure [J] shows the temporal decay of 
the density curve in a log-log plot, with a varied by increments of 10^'^ around 
0.62. The change of convexity is found between 0.627 and 0.629. The curves 
are obtained by averaging the data on Ng = 100 runs obtained on a ring size 
n = 20 000 and a samphng time of T = 200 000 steps. 

To improve our estimate of ac to a precision of 10~^, we extended the sam- 
pling time to T = 10^ and we increased the ring size to n = 40 000. Figured] 
shows the two values a~ = 0.6381 and a+ = 0.6383 for which we observed the 
change of convexity. We repeated the same type of progressive approximations 
of ttc for all the DP EGA. It was always possible to observe the change of con- 
vexity when varying a with an increment of 10~^. However, the sampling time 
T, the ring size n and the number of samples had to be adapted to each 
EGA to ensure a good stability of the measures (see below). The values of the 
estimated critical synchrony rates are reported in Table [1] 

4.1 Measurement of 5 

The value of 5 is given at the critical point by the slope of the density curve 
in a log- log scale (see Section 13. 4p . We report in Table [T] the time interval 
[Tmin, Tmax] uscd to computc 5. The lower limit Tmin is obtained by a rough 
estimate of the maximum transient time needed for the system to enter into 
the power law regime. The upper limit of the interval Tmax corresponds to the 
minimum time for which the deviation from a power law decay becomes visible 
for the curves a~ and that show the change of convexity. 

We bring to the reader's attention the fact that it is a difficult problem 
to estimate the error on 6. Indeed, besides the influence of noise, the value 
depends on the time interval [Tmin, Tmax] used to perform the fit. There is to 
our knowledge no general method for choosing this time interval. To estimate 
the error on (5, we computed the different values obtained when varying a to 
a~ and and by varying Tmin and Tmax by a factor 2. Numerical values of 

5 are reported in Table [1] for comparison with . Given our estimations of 
uncertainty, the results show good agreement with (J^p = 0.1595. 
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ECA 50 : Log-Log plot of d(t) for different values of alpha 
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ECA 50 : Log-Log plot of d(t) for different values of alpha 



f(t)=K.t'^-0.1595 
alpha= 0.6281 
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Figure 4: ECA 50: Determination of the critical synchrony rate ac-The straight 
line has slope —fop = —0.1595 and is plotted for reference. (above) Averages 
obtained on ~ 100 runs, ring size n = 20 000, sampling time T = 2 x 10'''. 
(below) Averages obtained on Ng = 50 runs and ring size n = 40 000, sampling 
time T = 10^. The arrow indicates the interval [Tmin, T^ax] used to perform the 
fit to measure S. 
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4.2 Stability of the measures 

An important question is to determine whether the visual estimations of the 
changes of convexity are satisfying and stable. A simple method to verify this 
is to estimate numerically the "linearity" of the density curves (in a log-log 
representation) . For the sake of simplicity, we used the root mean square error of 
the fits as a linearity estimator. For example, for EGA 50, fitting a power law in 
the time interval [300,30 000] for a" = 0.6181, = 0.6282, a+ = 0.6283 gives 
an error of 0.16, 0.10, 0.13 (dimensionless units), respectively, which confirms 
that the "most linear" curve is obtained with ac- 

The other point of care concerns the measurement of the critical exponent 5. 
Recall that we needed to distinguish the "power law" part [Tmin,Tmax] of the 
density curve and the "departure from power law" [Tniax,oo[, which can be 
convex or concave. To which extent does 5 depend on the choice of T^^^ and 
^max? Firstly, let us note that the length of transient time interval [0, Tmin] only 
depends on the EGA considered while the length of the "power law" part of the 
curve [Tmin, ?max] IS & function of the distance to the critical point |a — ad: the 
smaller this value, the longer the system follows the power law predictions. This 
means that the values of Tmax depend heavily on our choice of increment on a. 
In our experiment, we found out that it was necessary to take an increment on 
a as small as lO"** to ensure that the time interval [T,„in, Tmax] was large enough 
to measure 5 with a good precision. 

Secondly, recall that as we use finite size lattices, the system is metastable: 
it eventually reaches the absorbing state 0''', whatever the value of a. In order 
to infer the infinite-size limit from the simulations obtained on finite systems, 
we need to take n large enough to ensure that the trajectories do not that reach 
the absorbing state. For aU the DP EGA but EGA 18, we used n = 40 000 
to satisfy this condition ; for EGA 18, it was necessary to take n = 80 000 to 
prevent the system from reaching the absorbing state. 

Finally, the case of EGA 134 also needs to be underlined as its density curve 
has many inflexion points, which makes it harder to measure ckc and 5. 

5 Determination of (3 

In this second part of the experiment, we measure the (3 critical exponent by 
estimating the stationary density as a function of Aq, = |a — ad. 

5.1 Measurement of (3 

There are two different possibilities for varying Aq,: some authors use linear 
variation (e.g., |3]) while others use an exponential variation. As we expect 
the curve doo versus Aq, to be linear in a log- log plot, we varied Aq, with an 
exponential increment of \/2 from 8 x 10"'' to 0.512, which allowed us to obtain 
equally spaced points in the log-log scale. Figure [5] shows the evolution of the 
density versus time for different values of Aq . 



12 







alpha=0.6284 






alpha=0.6283 






alpha=0.6298 






alpha=0.6314 






alpha=0.6346 " 






alpha=0.6410 















Figure 5: EGA 50: Evolution of the density versus time for a > ac = 0.6282. 
The ring size is n = 2 x lO** ; averages use A^g = 10 runs with a samphng time 
T = 10^. 
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Figure 6: EGA 50: Determination of the; critical exponent P using the time decay 
properties (see text). The straight line shows theoretical prediction and has 
slope: /3dp = 0.2765. Note that both x and y axis are displayed in logarithmic 
scale. The arrow shows the fit interval used to calculate the /3 exponent. 
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Following the observations of the previous experiment, we estimated the 
stationary density doo{ct) by measuring d{t) during T = 10^ steps and by taking 
the average of d on the time interval [T/2,T]. To limit the influence of the 
noise, we repeated A's = 10 times the measure and took the average value. Note 
that even though T can be adjusted as a function of Aq, we prefer to take a 
fixed value for T for the sake of simplicity. This value was adjusted for the 
minimal value of , which makes it a priori suitable for larger values of Aq, 
(see Figure [5|). 

Figure [S] shows the estimated values of doo as a function of A^ . The visual 
comparison with the expected plot shows a good agreement with the predicted 
value of DP. For each EGA, we took t e [20 x 10"'', 20 x lO'^] as a fit interval. 
The estimated values of (3 are reported in Table [T]; as for the 6 exponent, they 
show good agreement with the DP prediction fS^p ^ 0.276. 

5.2 Stability of the measures 

To which extent do these measures depend on the setting of the parameters of 
the experiment? We observe on Figure [5] that the linear part of the density 
curve is limited by two competing phenomena. Note that the small values of 
Aq, are always overestimated. This can be explained by remarking that, near 
criticality the stationary density vanishes as: dao{a) ^ Aq,'^. Note that we also 
have d ~ near criticality, which implies that the time needed to get close to 
the stationary density increases exponentially with 1/Aq. This phenomenon, 
known as the critical slowing down {e.g., |17|). limits the measurement of the 
stationary density for small values of Aq,. 

For the higher values of Aq, the system "saturates" and no longer follows a 
power law (a density can not be higher than 1). The deviation from the power 
law is a phenomenon that is predicted by theory and that can be studied for 
its own interest. We prefer here to restrict our measures to the linear part of 
the curve. However, some authors advocate that the study of the deviations 
from the power law significantly improves the determination of a universality 
class [22]. To estimate the error on f3, we measured how its value changed when 
ac was varied in the interval [Q!~,a^] and when the fit interval was changed 
to [30 X 10"'', 30 X 10""^]. We noticed that the error was mainly due to the 
uncertainty on ac and had an order of magnitude of 10~^. 

6 The case of EGA 178 

For EGA 178, we repeated the protocol described above using the the kinks 
density instead of the density as an order parameter. Our best results were 
obtained with n = 80 000, T = 10* and iVg = 50: we located a change of 
convexity between a~ = 0.409 and = 0.411. However, by contrast with 
the previous DP EGA, this change of convexity was much more difficult to 
observe. This well-known difference is explained by the algebraic decay of the 
order parameter in PC or DP2 models. This produces straight lines in a log-log 
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Table 1: Numerical Results for the DP and DP2 rules : the digit between 
parentheses is uncertainty on the last digit. Compare with fop = 0.1595, 5dp2 = 
0.2856 and /3dp = 0.276. 



EGA 




5 


T ■ 

min 


T 

max 


f3 


6 


0.2825 (1) 


0.160 (3) 


1 X 10-^ 


3 X 10'^ 


0.27 (1) 


18 


0.71385 (5) 


0.157 (3) 


5 X 10^ 


2 X 10^ 


0.27 (1) 


26 


0.47485 (5) 


0.159 (1) 


1 X 10^ 


2 X 10^ 


0.27 (1) 


38 


0.04085 (5) 


0.160 (4) 


2 X lO'' 


4 X 10^ 


0.27 (1) 


50 


0.6282 (1) 


0.158 (3) 


3 X 10^ 


3 X 10* 


0.27 (1) 


58 


0.3398 (1) 


0.161 (4) 


1 X 10^ 


1 X 10^ 


0.28 (2) 


106 


0.8146 (1) 


0.157 (2) 


5 X 10^ 


3 X 10^ 


0.28 (1) 


134 


0.0821 (2) 


0.161 (5) 


2 X 10^ 


2 X 10*^ 


0.26 (3) 


146 


0.6751 (1) 


0.158 (3) 


5 X 10^ 


1 X 10^ 


0.27 (2) 


178 


0.410 (1) 


0.286 (4) 


1 X 10^ 


2 X 10^ 


?(?) 



scale which makes changes of convexity less visible than the exponential decay 
observed for the DP EGA. This explains why it is difficult to estimate ac with 
a better precision despite using a long computation time. 

In the time interval [100,20 000], the fit for the curve ac = 0.410 gives a 
slope of (5 = 0.2860, which is close to the predicted value (5dp2 ~ 0.2856. To 
estimate the error on 6, we repeated the fit with the same time interval, but 
with the curves obtained with a~ and q:+. We found that 5 varied less than 
4 X 10"'^. As the precision on 5 is good while the precision on ac is relatively 
poor for this rule, we did not measure the critical exponent /3. A possible way of 
improving our measure of the critical point ac consists in using a different exper- 
imental protocol, for example measuring the "dynamical" exponents obtained 
by starting from an initial condition close to the absorbing state. 

7 Discussion 

We experimentally investigated how qualitative changes of behaviour were trig- 
gered by gradual variations of the synchrony rate in asynchronous cellular au- 
tomata. The results show good evidence that these phenomena are second order 
phase transitions which belong to the directed percolation (DP) and parity con- 
servation (PC- DP2) universality classes. 

From a practical point of view, the main limitation for identifying the phase 
transitions was the huge amount of computation time required for measuring 
the critical exponents : for each EGA, we used approximately lO'^^ computa- 
tions of the local rule in order to determine these exponents with two or three 
digits. This task represented several months of calculus with several personal 
computers. It calls for further studies on how to improve the compiitation pro- 
cess {e.g., the generation of random numbers) or on how to take advantage of 
the use of massively parallel computing devices. 
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7.1 The origin of phase transitions 

The observation of the synchronous behaviour of the rules we studied indicate 
that there is certainly no straightforward relation with the existing classifica- 
tions. For example, EGA 50 is "periodic" (or Wolfram class II) while EGA 18 
is "chaotic" (or Wolfram class III). This indicates that some particular asyn- 
chronous cellular models obey the same macroscopic laws near their transition 
points despite of having different definitions at the microscopic cell-scale. In 
other words, the asynchronous updating of cellular systems unveils another type 
of complexity which still needs to be understood. 

A first path for discussing these results is to consider the famous conjec- 
ture by Janssen and Grassberger {e.g., see [17] for a short presentation). It 
states that a model should belong to the DP universality class if it satisfies the 
following criteria: (a) the model displays a continuous phase transition from 
a fluctuating active phase to a unique absorbing state, (b) it is possible to 
characterise the phase transition by a positive one-component order parameter, 
(c) the dynamics of the model is defined by short-range process, and (d) there 
exists no additional symmetries or quenched randomness (i.e., small and stable 
topological modifications) . 

For all the DP EGA, we saw that condition (a) was fulfilled with 0-^ as the 
absorbing state. It is interesting to note that configuration O'^ is a fixed point 
for all the DP rules, but EGA 134 and 146 also have I''' as a fixed point. This 
confirms that the informal notion of "absorbing state" can not be trivially iden- 
tified with the fixed point mathematical property. Gondition (b) was fulfilled 
by the density or by the kinks density. Gondition (c) is true by definition of 
cellular automata. 

Analysing condition (d) is more interesting since we see a difference between 
space symmetry, which is possessed by EGA 50 for example, and state symmetry 
{i.e., invariance under and 1 exchanging), which is absent for all the DP 
EGA. This may also explain why EGA 178, which has both symmetries, has a 
peculiar behaviour and was found in the PC-DP2 universality class. Note that 
conditions (b), (c), (d) can be easily verified but a challenging question consists 
in explaining why only a small fraction of the 256 EGA also verify condition (a) 
and thus exhibit DP behaviour. 

7.2 Perspectives 

At this point, it is still an intriguing question to understand the origin of out- 
of-equilibrium phase transitions in cellular automata, and more generally, in 
particle systems. The models we exhibited count among the simplest models 
that display out-of-equilibrium phase transitions. In that respect, they can be 
used as a good basis for exploring the origin of phase transitions by analytical 
means or by numerical simulations. 

A further step for understanding asynchronous cellular automata is to make 
a "reduction" between the DP rules, i.e., to show that if one of them exhibits 
such a phase transition, then the others behave similarly. Another step is to ex- 
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tend such a reduction to other well-studied systems such as the Domany-Kinzel 
probabilistic C A [8] . This implies uniting the two types of phase transitions we 
observed: recall that going from the inactive to active phase was obtained either 
by increasing (DPhi EGA) or by decreasing (DPiow EGA) the control parameter. 
To our knowledge, it is the first example where these two-way transitions are 
observed. 

A more ambitious possibility for unifying the study of all the DP EGA is to 
consider the larger set of probabilistic elementary cellular automata, which in- 
clude the a-asynchronous EGA. In this space, which is homeomorphic to [0, 1]®, 
the problem is to determine whether the phase transitions can be explained in 
terms of crossing of a hypersurface (see also |7])- 

Our view is that there exists a wide range of problems that could benefit 
from the study of cellular systems with out-of-equilibrium phase transitions. For 
example, in biology, can we explain the trigger of the self-organisation phase in 
cellular societies by using similar models [T2l[2|? In an engineering context, can 
we use phase transitions to design systems that change their behaviour without 
any external control ? 
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